package org.korifey.ffalg;

/**
 * Number Theory.
 * 
 * @author daivanov
 */
public final class NumberTheory {

	/**
	 *
	 */
	private NumberTheory() {
	}

	/**
	 * 
	 * @param x
	 * @return number of significant bits in x
	 */
	public static int nbits(long x) {
		assert x >= 0;

		int res = 0;
		while (x != 0) {
			res += 1;
			res >>= 1;
		}
		return res;
	}

	/**
	 * 
	 * @param x
	 * @return
	 */
	public final static int upperLimitOfPrimes(long x) {
		if (x <= 10) {
			return 4;
		} else {
			double p = x / Math.log(x);
			double res;
			if (x <= 1000) {
				res = p * 1.3;
			} else if (x <= 10000) {
				res = p * 1.183;
			} else if (x <= 100000) {
				res = p * 1.136;
			} else if (x <= 1000000) {
				res = p * 1.11;
			} else if (x <= 10000000) {
				res = p * 1.085;
			} else if (x <= 100000000) {
				res = p * 1.0713;
			} else if (x <= 1000000000) {
				res = p * 1.0614;
			} else {
				res = p * 1.0538;
			}

			return (int) res;
		}
	}

	/**
	 * 
	 * @param x
	 * @return
	 */
	public final static int probableNumberOfPrimes(long x) {
		assert x > 0;
		return (int) (x / Math.log(x));
	}

	/**
	 * 
	 */
	private static int opcount = 0;

	/**
	 * 
	 * @param max
	 * @return
	 */
	public final static int[] erathosphen(int max) {
		assert max > 0;
		int n = 0;

		FFBitSet bitSet = new FFBitSet(max + 1); // don't count 0 and 1
		int i;
		opcount = 0;

		for (i = 2; i <= (int) Math.sqrt(max); i++) {
			// opcount++;
			if (!bitSet.get(i)) {
				n++;
				for (int j = i * i; j <= max; j += i) {
					bitSet.set(j);
					// opcount ++;
				}
			}
		}
		for (; i <= max; i++) {
			// opcount++;
			if (!bitSet.get(i)) {
				n++;
			}
		}

		int[] res = new int[n];
		int cur = 0;
		for (i = 2; i <= max; i++) {
			if (!bitSet.get(i)) {
				res[cur++] = i;
			}
		}

		return res;
	}

	public final static boolean isPrime(int x, int[] eratosphen) {
		assert x >= 0;
		if (x <= 1)
			return false;
		for (int i = 0; i < eratosphen.length; i++) {
			if (x % eratosphen[i] == 0)
				return false;
		}
		return true;
	}

	/**
	 * 
	 * @param a
	 * @param b
	 * @return
	 */
	private final static long gcdOfTwo(long a, long b) {
		assert a != 0;
		assert b != 0;
		a = Math.abs(a);
		b = Math.abs(b);

		long c = a % b;
		while (c != 0) {
			a = b;
			b = c;
			c = a % b;
		}
		return b;
	}

	/**
	 * 
	 * @param a
	 * @param b
	 * @return
	 */
	public final static long gcd(long first, long... rest) {
		long res = first;
		for (long x : rest) {
			res = gcdOfTwo(res, x);
		}
		return res;
	}

	public final static long gcd(long[] arg) {
		assert arg.length > 0;
		long res = arg[0];
		return gcd(res, arg);
	}

	/**
	 * 
	 * @param a
	 * @param b
	 * @return
	 */
	public final static long gcm(long a, long b) {
		return (a * b) / gcd(a, b);
	}

	/**
	 * 
	 * @author daivanov
	 * 
	 */
	public static class GcdExtResult {
		private long x;
		private long y;
		private long gcd;

		/**
		 * 
		 * @param x
		 * @param y
		 * @param gcd
		 */
		public GcdExtResult(long x, long y, long gcd) {
			super();
			this.x = x;
			this.y = y;
			this.gcd = gcd;
		}

		/**
		 * 
		 * @return
		 */
		public long x() {
			return x;
		}

		/**
		 * 
		 * @return
		 */
		public long y() {
			return y;
		}

		/**
		 * 
		 * @return
		 */
		public long gcd() {
			return gcd;
		}

		@Override
		public int hashCode() {
			final int prime = 31;
			int result = 1;
			result = prime * result + (int) (gcd ^ (gcd >>> 32));
			result = prime * result + (int) (x ^ (x >>> 32));
			result = prime * result + (int) (y ^ (y >>> 32));
			return result;
		}

		@Override
		public boolean equals(Object obj) {
			if (this == obj)
				return true;
			if (obj == null)
				return false;
			if (getClass() != obj.getClass())
				return false;
			GcdExtResult other = (GcdExtResult) obj;
			if (gcd != other.gcd)
				return false;
			if (x != other.x)
				return false;
			if (y != other.y)
				return false;
			return true;
		}

		@Override
		public String toString() {
			return "GcdExtResult [x=" + x + ", y=" + y + ", gcd=" + gcd + "]";
		}

	}

	public final static GcdExtResult gcdExt(long a, long b) {
		assert a > 0;
		assert b > 0;

		long x0 = 1, x1 = 0, r0 = a, r1 = b;

		while (true) {
			long r2 = r0 % r1;
			if (r2 == 0) {
				break;
			}

			// we can compute only first coefficient, second one we will
			// calculate at the and
			long q = -r0 / r1;
			long x2 = q * x1 + x0;

			// move to the next iteration
			r0 = r1;
			r1 = r2;
			x0 = x1;
			x1 = x2;
		}
		return new GcdExtResult(x1, (r1 - x1 * a) / b, r1);
	}

	/**
	 * ax=b(mod c)
	 * 
	 * @param a
	 * @param b
	 * @param m
	 */
	public final static AffineFunction solveModEquation(long a, long b, long c) {
		assert a >= 0;
		assert b >= 0;
		assert c > 0;

		a = a % c;
		b = b % c;

		if (a == 0) {
			if (b == 0) {
				return AffineFunction.linearFunction(1);
			} else {
				return null;
			}
		} else {
			if (b == 0) {
				return AffineFunction.linearFunction(gcm(a, c) / a);
			} else {
				GcdExtResult r = gcdExt(a, c);
				if (b % r.gcd() != 0) {
					return null;
				} else {
					GcdExtResult gcdRes = gcdExt(a, c);
					long offset = gcdRes.x * (b / gcdRes.gcd);
					offset = (offset % c + c) % c;
					return AffineFunction.newInstance(offset, c / gcdRes.gcd);
				}
			}
		}
	}

	public final static long divMod(long a, long b, long mod) {
		AffineFunction f = solveModEquation(b, a, mod);
		if (f == null)
			return 0;
		else {
			long coeff = (a - b * f.getOffset()) / (b * f.getCoeff()[0]);
			return (f.getOffset() + coeff * f.getCoeff()[0]) % mod;
		}

	}

	public final static long inverseMod(long a, long mod) {
		assert a > 0;
		assert mod > 0;

		AffineFunction modRes = solveModEquation(a, 1, mod);
		return modRes == null ? 0 : modRes.getOffset();
	}

	public static long powMod(long x, long p, long mod) {
		assert x > 0;
		assert p > 0; // TODO
		assert mod > 0;

		long res = 1;
		x = x % mod;

		long factor = x;
		while (p > 0) {
			if ((p & 1) != 0)
				res = (res* factor) % mod;
			factor = (factor*factor) % mod;
			p >>= 1;
		}
		return res;
	}
	
	public class FactorizeResult {		
		int[] primeNumbers;
		int[] primeCounts;
		int lastFactor;
	}
	
	public static FactorizeResult factorize(int x, int[] primes) {
		assert (x > 0);
		int[] factors = new int[64];
		int cur = 0;		
		
		for (int i=0; primes[i]*primes[i] <= x; i++) {
			while (x % primes[i] == 0) {
				x/=primes[i];
				factors[cur++] = i;
			}
		}
		return null;
	}

	public static void main(String[] args) {

		testDist();
		if (true)
			return;

		int i = 100010000;
		int prob = probableNumberOfPrimes(i);
		long start = System.currentTimeMillis();
		int[] primes = erathosphen(i);
		double act = primes.length;
		System.out
				.println(String
						.format("N = %d : primes(n) = %d, DIFF = %.4f : TIME = %d : opcount = %d",
								i, primes.length, (act / prob),
								(System.currentTimeMillis() - start), opcount));
		int cur = 0;
		for (int j = 2; j < i; j++) {
			if (cur < primes.length && primes[cur] == j) {
				cur++;
				prob = upperLimitOfPrimes(j);
				if (prob < cur) {
					System.out.println("Bad prime: " + j + " prob: " + prob
							+ ", act: " + cur + ", ratio: "
							+ (cur / (double) probableNumberOfPrimes(j)));
				}
			}
		}
	}

	public static void testDist() {
		int[] dists = new int[100];
		int[] primes = erathosphen(1000000000);
		
		int longdists = 0;
		long e = 0;
		for (int i=1; i<primes.length; i++) {
			int dst = primes[i] - primes[i-1];
			e+=dst;
			if (dst >= dists.length) {
				longdists++;
				continue;
			}			
			dists[dst] ++;
		}
		
		System.out.println("Average dist =" + e / (primes.length-1));
		for (int i=0; i<dists.length; i++) {
			if (dists[i] != 0) System.out.println(i+":"+dists[i]);			
		}
		System.out.println("long:"+longdists);
	}

}
